! This module contains coefficients

module coeffpg
   use commirr

   integer, parameter:: ne  = 169  &
                      , nee = 201  &
					  ,	np  = 421

   real(8), dimension(nee):: xxee, yyee, ypee	           ! pair production

   real(8), dimension(ne ):: rr  , B   , D   , BP, DP, DPP !

   real(8), dimension(np ):: dlgE					 &	   ! lg(E/eV)
                           , Ppp , Ppn , Pnn , Pnp   &     ! Probabilities
                           , Bpp , Bpn , Bnn , Bnp   &     ! Relat. energy losses 
                           , DBpp, DBpn, DBnn, DBnp        ! Deriv. rel.ener.losses


   save
     
!--------------------------------------------------------
                    contains
!--------------------------------------------------------

!***************************************************
! Read coefficients for e^+ + e^- and              *
! nucleon photoproduction in CMBR  at z=0          *
!***************************************************
   Subroutine ReadCoeff
      implicit none
	  integer:: k
!-------------
      open(1,file='Data/BBPee.dat',status='old')
	     read(1,*)
	     read(1,*)
         read(1,'(f8.4,f18.11,e20.11)')   &
            (xxee(k), yyee(k), ypee(k),k=1,nee)
	  close(1)
!------------- 
      open(1,file='Data/B5.dat',status='old')
	     read(1,*)
         read(1,'(F8.4,2F18.11,3E20.11)') &
            (rr(k),B(k),D(k),BP(k),DP(k),DPP(k),k=1,ne)
      close(1)
!-------------
      open(1,file='Data/ProbNN.dat',status='old') 
         read(1,'(5/,(5F13.7))')   &
            (dlgE(k),Ppp(k),Ppn(k),Pnn(k),Pnp(k),k=1,np)
      close(1)    
!-------------
      open(1,file='Data/BetaNN.dat',status='old') 
         read(1,'(5/,(5F13.7))')  &
            (dlgE(k),Bpp(k),Bpn(k),Bnn(k),Bnp(k),k=1,np)
      close(1)    
!-------------
      open(1,file='Data/DrvBeta.dat',status='old') 
         read(1,'(5/,(F13.7,4E15.7))')	&
        (dlgE(k),DBpp(k),DBpn(k),DBnn(k),DBnp(k),k=1,np)
      close(1)    
!-------------
      write(*,*) ' Data is read !'

   end Subroutine ReadCoeff
!***************************************************


!***************************************************
! Calculates all approximations to the             *
! coefficients  of e^+ + e^- and                   *
! nucleon photoproduction in CMBR  at z=0          *
!***************************************************
!  r = dlog10(E) in eV                             *
!***************************************************
   Subroutine AprBD(r,Bee,BPee,Bout,Dout,BPout,DPout,DPPout)
	  use msimsl, only: DCSIEZ
	  use prmmirr, only: Flin

      implicit none

	  real(8), dimension(1), intent( in):: r
	  real(8), dimension(1), intent(out):: Bee, BPee, Bout  &
	                                     , Dout, BPout, DPout, DPPout
      real(8):: rv
	  rv = r(1)

!**********************************************************************
      if(rv.le.xxee(1)) then
         Bee(1)  = flin(rv,xxee(1),xxee(2),yyee(1),yyee(2))   ! lg B^ee
         BPee(1) = flin(rv,xxee(1),xxee(2),ypee(1),ypee(2))   ! dB^ee/dlnE 
      elseif(rv.gt.xxee(nee)) then
         Bee(1)  = flin(rv,xxee(nee-1),xxee(nee),yyee(nee-1),yyee(nee))
         BPee(1) = flin(rv,xxee(nee-1),xxee(nee),ypee(nee-1),ypee(nee))
      else
         CALL DCSIEZ (nee, xxee, yyee,1,r,Bee )
         CALL DCSIEZ (nee, xxee, ypee,1,r,BPee)
      endif
!**********************************************************************
      if(rv.le.rr(1)) then
         Bout(1)   = flin(rv,rr(1),rr(2),B  (1),B  (2))
         Dout(1)   = flin(rv,rr(1),rr(2),D  (1),D  (2))
         BPout(1)  = flin(rv,rr(1),rr(2),BP (1),BP (2))
         DPout(1)  = flin(rv,rr(1),rr(2),DP (1),DP (2))
         DPPout(1) = flin(rv,rr(1),rr(2),DPP(1),DPP(2))
      elseif(rv.gt.rr(ne)) then
         Bout(1)   = flin(rv,rr(ne-1),rr(ne),B  (ne-1),B  (ne))
         Dout(1)   = flin(rv,rr(ne-1),rr(ne),D  (ne-1),D  (ne))
         BPout(1)  = flin(rv,rr(ne-1),rr(ne),BP (ne-1),BP (ne))
         DPout(1)  = flin(rv,rr(ne-1),rr(ne),DP (ne-1),DP (ne))
         DPPout(1) = flin(rv,rr(ne-1),rr(ne),DPP(ne-1),DPP(ne))
      else
         CALL DCSIEZ (ne, rr, B  ,1,r,Bout  )
         CALL DCSIEZ (ne, rr, D  ,1,r,Dout  )
         CALL DCSIEZ (ne, rr, BP ,1,r,BPout )
         CALL DCSIEZ (ne, rr, DP ,1,r,DPout )
         CALL DCSIEZ (ne, rr, DPP,1,r,DPPout)
      endif
!**********************************************************************
!      Bout = 1d1**Bout
!      Dout = 1d1**Dout
!**********************************************************************
   end Subroutine AprBD


!**********************************************************************
   Subroutine AprCoefZ(z1,dlnE)
	  use prmmirr, only: Lambda, Omega_m, ln10
	  use commirr, only: Xmir, C
       
	  implicit none

	  integer:: k
 	  real(8), intent(in):: z1, dlnE
	  real(8):: z13, fact_z
	  real(8), dimension(1):: r, Bee, BPee,Bout  &
	                        , Dout, BPout, DPout, DPPout

	  z13    = z1*z1*z1
	  fact_z = z13/dsqrt(Omega_m*z13 + Lambda)
	  r(1)   = dlog10(z1) + dlnE/ln10   ! lg((1+z)*E)      

	  call AprBD(r,Bee,BPee,Bout,Dout,BPout,DPout,DPPout)
	  
	  C(1) = 1d1**(Bee(1) + dlog10(1d0 + 1d1**(Bout(1) - Bee(1))))
	  C(2) = 1d1**Dout(1)
	  C(3) = BPee(1) + BPout(1)
	  C(4) = DPout(1)
	  C(5) = DPPout(1)

	  do k=1,5
	     C(k) = fact_z*C(k)
	  enddo	
	  
   end Subroutine AprCoefZ
!**********************************************************************
      

!**********************************************************************
   Subroutine AprCoefZMir(z1,dlnE)
	  use commirr, only: Xmir, Xmir3, Cmir, Lambda, Omega_m, ln10
       
	  implicit none

	  integer:: k
	  real(8), intent(in):: z1, dlnE
	  real(8):: z13, fact_z, Emir
	  real(8), dimension(1):: r, Bee, BPee,Bout  &
	                        , Dout, BPout, DPout, DPPout

	  z13 = z1*z1*z1
	  fact_z = z13/dsqrt(Omega_m*z13+Lambda)

      Emir = dexp(dlnE) 	
	  
	  r(1)    = dlog10(z1) + (dlnE + dlog(Xmir))/ln10   ! lg((1+z)*E)      
  
	  call AprBD(r,Bee,BPee,Bout,Dout,BPout,DPout,DPPout)

	  Cmir(1) = 1d1**(Bee(1) + dlog10(1d0 + 1d1**(Bout(1) - Bee(1))))
	  Cmir(2) = 1d1**Dout(1)
	  Cmir(3) = BPee(1) + BPout(1)
	  Cmir(4) = DPout(1)
	  Cmir(5) = DPPout(1)

      do k=1,5
	     Cmir(k) = fact_z*Cmir(k)*Xmir3
	  enddo	
	  
   end Subroutine AprCoefZMir
!**********************************************************************

end module coeffpg